*DECK ARCLEN
      SUBROUTINE ARCLEN(XMIN,XMAX,A,B,XN,IMAX,S)
      DX=(XMAX-XMIN)/IMAX
      S=0.0
      X=XMIN-DX
      DO 500 I=1,2100
      XS=X
      X=X+DX
      IF(X.GT.XMAX) DX=XMAX-XS
      IF(X.GT.XMAX) X=XMAX
      YS=Y
      CALL FCT(A,B,XN,1,X,Y)
      IF(I.EQ.1) GO TO 500
      D=SQRT(DX**2+(Y-YS)**2)
      S=S+D
      IF(X.GE.XMAX) RETURN
  500 CONTINUE
      WRITE(6,600) IMAX,XN,DX,S,D,X,XS,Y,YS,XMIN,XMAX,A,B
  600 FORMAT(29H STOPPED IN SUBROUTINE ARCLEN /I5,12E10.5)
      STOP
      END
*DECK AXIS
      SUBROUTINE AXIS(XMIN,XMAX,YMIN,YMAX,NXINTS,NYINTS,XRN,XRX,YRN,YRX)
C
C        THIS SUBROUTINE IS ONE OF THOSE THAT PLOTS GRIDS
C        ON THE SC4020 PLOTTER.
C
C     THIS SUBROUTINE DRAWS THE X AND Y AXES, MARKS OFF THE
C     INTERVALS, SETS THE SCALE VALUES, AND MOVES THE FILM TO A NEW
C     FRAME.
C
C     XMIN IS THE MINIMUM VALUE OF THE X (HORIZONTAL) VARIABLE.
C     XMAX IS THE MAXIMUM VALUE OF THE X VARIABLE.
C     YMIN IS THE MINIMUM VALUE OF THE Y (VERTICAL) VARIABLE.
C     YMAX IS THE MAXIMUM VALUE OF THE Y VARIABLE.
C     NXINTS IS THE NUMBER OF INTERVALS INTO WHICH THE X AXIS IS TO
C          BE DIVIDED.
C     NYINTS IS THE NUMBER OF INTERVALS INTO WHICH THE Y AXIS IS TO
C          BE DIVIDED.
C
C     DO INITIALIZATION
      CALL SMALLV
      CALL FRAMEV(3)
      NX=NXINTS
      NY=NYINTS
C     IF NX OR NY .LE. 0 MAKE THEM EQUAL TO 1
      IF (NX .LE. 0) NX=1
      IF (NY .LE. 0) NY=1
C     SET UP VALUES TO BE USED FOR MARGINS
      ML=123
      MR=923
      MB=123
      MT=923
C     DRAW X AXIS
      CALL LINEV(ML,MB,MR,MB)
C     DRAW Y AXIS
      CALL LINEV(ML,MB,ML,MT)
C     DETERMINE INCREMENTS FOR TIC MARKS
      DX=(XMAX-XMIN)/NX
      DY=(YMAX-YMIN)/NY
C     SCALE X AND Y VALUES
      CALL XSCALV(XMIN,XMAX,ML,100)
      CALL YSCALV(YMIN,YMAX,MB,100)
C     DRAW TIC MARKS ON THE X AXIS
      DX=(XRX-XRN)/NX
      CALL LINRV(1,MB-20,MB-5,MB+5,XMIN,XMAX,DX,0,-1,3,8)
C     CALL LINEV(NXV(SING),NYV(0.0)-20,NXV(SING),NYV(0.0))
C     DRAW TIC MARKS ON Y AXIS
      DY=(YRX-YRN)/NY
      CALL LINRV(2,ML- 90,ML-5,ML+5,YRN,YRX,DY,0,-1,3,10)
      RETURN
      END
*DECK BDRYS
      SUBROUTINE BDRYS(JMIN,JMAX,KMIN,KMAX,XMIN,XMAX,YMIN,YMAX,XN,DX1,
     1                   DY1,DS1,A,B,X,Y,IB,D1,D2,WL,WR,WB,WT,WOER,WOET)
      DIMENSION X(40,40),Y(40,40),WL(1),WR(1),WB(1),WT(1),WOER(1),
     1                                                    WOET(1),IB(10)
      EPSLON=.001
      J1=JMIN+1
      K1=KMIN+1
      KMAXM1=KMAX-1
      JMAXM1=JMAX-1
C     LEFT SIDE OF RECTANGULAR GRID -- IB(5 OR 6)
      IBI=5
      IF(KMIN.GT.1) IBI=6
      IF(IB(IBI).EQ.1) CALL WCALC(X,Y,WL,1,1,K1,KMAXM1)
      IF(IB(IBI).EQ.1) GO TO 300
      X(1,KMIN)=0.0
      Y(1,KMIN)=YMIN
      X(1,KMAX)=0.0
      Y(1,KMAX)=YMAX
      EY=EPSIL(YMAX-YMIN,0.0,DY1,KMAX-KMIN+1,0.001,100,1,1)
      IF(KMIN.GT.1) GO TO 230
      WLK=(KMAX-KMIN-1)*ALOG(1.0+EY)/(KMAX-KMIN)
      DO 200 KK=K1,KMAXM1
      K=KMAXM1-KK+K1
      X(1,K)=X(1,K+1)
      Y(1,K)=Y(1,K+1)-DY1*(1.0+EY)**(KMAXM1-K)
      WL(K)=WLK
  200 CONTINUE
      GO TO 300
  230 DO 240 K=1,KMIN
  240 WL(K)=0.0
      EYK=1.0/(1.0+EY)-1.0
      WLK=(KMAX-KMIN-1)*ALOG(1.0+EYK)/(KMAX-KMIN)
      DO 250 K=K1,KMAXM1
      X(1,K)=X(1,K-1)
      Y(1,K)=Y(1,K-1)+DY1*(1.0+EY)**(K-K1)
      WL(K)=WLK
  250 CONTINUE
  300 IF(IB(1).NE.4) GO TO 350
C     BOTTOM, RIGHT SIDE AND TOP OF RECTANGULAR GRID BY ELLIPSE TEST CASE
C
      CALL ELIPSE(XMAX,YMAX,X,Y,JMAX,KMAX,IB,DX1,WR,WB,WT,A,B,XN,EPSLON)
      GO TO 700
C     BOTTOM OF RECTANGULAR GRID -- IB(1 OR 2)
  350 IBI=1
      IF(JMIN.GT.1) IBI=2
      IF(IB(IBI).EQ.1) CALL WCALC(X,Y,WB,J1,JMAXM1,1,1)
      IF(IB(IBI).EQ.1) GO TO 480
      E=EPSIL(XMAX-XMIN,0.0,DX1,JMAX-JMIN+1,0.001,100,1,1)
      X(JMAX,1)=XMAX
      Y(JMAX,1)=0.0
      X(JMIN,1)=XMIN
      Y(JMIN,1)=0.0
      IF(JMIN.GT.1) GO TO 430
      WBJ=(JMAX-JMIN-1)*ALOG(1.0+E)/(JMAX-JMIN)
      DO 400 JJ=J1,JMAXM1
      J=JMAXM1-JJ+J1
      X(J,1)=X(J+1,1)-DX1*(1.0+E)**(JMAXM1-J)
      Y(J,1)=Y(J+1,1)
      WB(J)=WBJ
  400 CONTINUE
      GO TO 480
  430 DO 440 J=1,JMIN
  440 WB(J)=0.0
      EJ=1.0/(1.0+E)-1.0
      WBJ=(JMAX-JMIN-1)*ALOG(1.0+EJ)/(JMAX-JMIN)
      DO 450 J=J1,JMAXM1
      X(J,1)=X(J-1,1)+DX1*(1.0+E)**(J-J1)
      Y(J,1)=Y(J-1,1)
      WB(J)=WBJ
  450 CONTINUE
C     TOP AND RIGHT SIDE OF SUPER ELLIPSE -- IB(7&8 OR 9&10)
  480 JMX=JMAX
      IF(JMIN.GT.1) JMX=JMIN
      KMX=KMAX
      IF(KMIN.GT.1) KMX=KMIN
      IB7=7
      IF(KMIN.EQ.1) IB7=9
      IB8=8
      IF(JMIN.EQ.1) IB8=10
      IF(IB7.EQ. 7.AND.IB(IB7).EQ.1)CALL WCALC(X,Y,WOER,JMX,JMX,2,KMX-1)
      IF(IB7.EQ. 9.AND.IB(IB7).EQ.1)CALL WCALC(X,Y,WR  ,JMX,JMX,2,KMX-1)
      IF(IB8.EQ. 8.AND.IB(IB8).EQ.1)CALL WCALC(X,Y,WOET,2,JMX-1,KMX,KMX)
      IF(IB8.EQ.10.AND.IB(IB8).EQ.1)CALL WCALC(X,Y,WT  ,2,JMX-1,KMX,KMX)
      IF(IB(IB7).EQ.1.AND.IB(IB8).EQ.1) GO TO 700
      CALL SLOPE(A,B,XN,XD,YD)
      CALL ARCLEN(0.0,XD,A,B,XN,512,S1)
      CALL ARCLEN(0.0,YD,B,A,XN,512,S2)
      E1=EPSIL(S1,0.0,DS1,JMX,0.001,100,1,1)
      E2=EPSIL(S2,0.0,DS1,KMX,0.001,100,1,1)
      WTJ=(JMX-2)*ALOG(1.0+E1)/(JMX-1.0)
      WRK=(KMX-2)*ALOG(1.0+E2)/(KMX-1.0)
      DX=DS1*(1.0+E1)**(JMX-2)
      IF(IB(IB8).EQ.1) GO TO 550
      X1=0.0
      DO 500 J=2,JMX
      DS=DS1*(1.0+E1)**(JMX-J)
      CALL XOFS(DS,X1,DX,A,B,XN,EPSLON)
      X(J,KMX)=X1
      CALL FCT(A,B,XN,1,X1,Y(J,KMX))
      WT(J)=WTJ
      IF(JMIN.GT.1) WOET(J)=WTJ
  500 CONTINUE
  550 IF(IB(IB7).EQ.1) GO TO 700
      Y1=0.0
      DY=DX
      DO 600 K=2,KMX
      DS=DS1*(1.0+E2)**(KMX-K)
      CALL XOFS(DS,Y1,DY,B,A,XN,EPSLON)
      CALL FCT(A,B,XN,2,X(JMX,K),Y1)
      Y(JMX,K)=Y1
      WR(K)=WRK
      IF(KMIN.GT.1) WOER(K)=WRK
  600 CONTINUE
  700 IF(JMIN.GT.1) GO TO 800
      D1=Y(1,2)-Y(1,1)
      D2=Y(JMAX,2)-Y(JMAX,1)
      CALL WCALC(X,Y,WR,JMAX,JMAX,2,KMAXM1)
      CALL WCALC(X,Y,WL,1,1,2,KMAXM1)
      CALL WCALC(X,Y,WT,2,JMAXM1,KMAX,KMAX)
      CALL WCALC(X,Y,WB,2,JMAXM1,1,1)
      RETURN
C     RIGHT OUTER BOUNDARY -- IB(3)
  800 IF(IB(3).EQ.1) CALL WCALC(X,Y,WR,JMAX,JMAX,2,KMAXM1)
      IF(IB(3).EQ.2) CALL ROB(X,Y,WR,JMIN,JMAX,KMIN,KMAX,YMAX)
      IF(IB(3).GT.0) GO TO 1000
      DY=YMAX/KMAXM1
      DO 900 K=2,KMAX
      X(JMAX,K)=X(JMAX,K-1)
      Y(JMAX,K)=Y(JMAX,K-1)+DY
  900 WR(K)=0.0
C     TOP OUTER BOUNDARY -- IB(4)
 1000 IF(IB(4).EQ.1) CALL WCALC(X,Y,WT,2,JMAXM1,KMAX,KMAX)
      IF(IB(4).EQ.2) CALL TOB(X,Y,WT,JMIN,JMAX,KMIN,KMAX)
      IF(IB(4).GT.0) RETURN
      DX=XMAX/JMAXM1
      DO 1100 J=2,JMAX
      X(J,KMAX)=X(J-1,KMAX)+DX
      Y(J,KMAX)=Y(J-1,KMAX)
 1100 WT(J)=0.0
      RETURN
      END
      SUBROUTINE ELIPSE(XMAX,YMAX,X,Y,JMAX,KMAX,IB,DX1,WR,WB,WT,A,B,XN,
     !                                                           EPSLON)
      COMMON /PLOTC/ SING
C     THIS SUBROUTINE SETS UP BOTTOM AND TOP OF ELLIPSE TEST CASE
C     THAT IS IB(1),IB(9),IB(10)
      DIMENSION X(40,40),Y(40,40),IB(1),WR(1),WB(1),WT(1)
      SING=SQRT(XMAX**XN-YMAX**XN)
      KMAXM1=KMAX-1
      JMAXM1=JMAX-1
      X(JMAX,1)=SING
      Y(JMAX,1)=0.0
      X(1,1)=0.0
      Y(1,1)=0.0
      X(JMAX,KMAX)=XMAX
      Y(JMAX,KMAX)=0.0
      DX=DX1
      E=EPSIL(XMAX-SING,0.0,DX,KMAX,0.001,100,1,1)
      WRK=(KMAX-2)*ALOG(1.0+E)/KMAXM1
C     NEXT 3 STATEMENTS CAUSE EQUAL SPACING
C     E=0.0
C     DX=(XMAX-SING)/KMAXM1
C     WRK=0.0
      DO 100 KK=2,KMAXM1
      K=KMAXM1-KK+2
      X(JMAX,K)=X(JMAX,K+1)-DX*(1.0+E)**(KMAXM1-K)
      Y(JMAX,K)=0.0
      WR(K)=WRK
  100 CONTINUE
      CALL SLOPE(A,B,XN,XD,YD)
      CALL ARCLEN(0.0,XD,A,B,XN,512,S1)
      CALL ARCLEN(0.0,YD,B,A,XN,512,S2)
      DS1=.5*(S1+S2)/JMAXM1
      E=EPSIL(S1+S2,0.0,DS1,JMAX,0.001,100,1,1)
      WTJ=(JMAX-JJ-1)*ALOG(1.0+E)/(JMAX-JJ)
      DX=DS1
      X1=0.0
C     NEXT 3 STATEMENTS CAUSE EQUAL SPACING
C     E=0.0
C     DS1=2.0*DS1
C     WTJ=0.0
      DO 300 J=2,JMAXM1
      DS=DS1*(1.0+E)**(JMAX-J)
      CALL XOFS(DS,X1,DX,A,B,XN,EPSLON)
      IF(X1.GT.XD) GO TO 400
      X(J,KMAX)=X1
      CALL FCT(A,B,XN,1,X1,Y(J,KMAX))
      WT(J)=WTJ
  300 CONTINUE
      RETURN 0
  400 JJ=J
      Y1=0.0
      DY=DS1
      DO 500 J=JMAXM1,JJ,-1
      DS=DS1*(1.0+E)**(JMAXM1-J)
      CALL XOFS(DS,Y1,DY,B,A,XN,EPSLON)
      Y(J,KMAX)=Y1
      CALL FCT(A,B,XN,2,X(J,KMAX),Y1)
      WT(J)=WTJ
  500 CONTINUE
      DXX=X(2,KMAX)-X(1,KMAX)
      E=EPSIL(SING,0.0,DXX,JMAX,0.001,100,1,1)
      WBJ=(JMAX-2)*ALOG(1.0+E)/JMAXM1
C     NEXT 3 STATEMENTS CAUSE EQUAL SPACING
C     E=0.0
C     DXX=SING/JMAXM1
C     WBJ=0.0
      DO 600 J=2,JMAXM1
      X(J,1)=X(J-1,1)+DXX*(1.0+E)**(J-2)
      Y(J,1)=0.0
      WB(J)=WBJ
  600 CONTINUE
      RETURN
      END
*DECK EPSIL
      FUNCTION EPSIL(FMX,FMIN,DFM,NPT,FPCC,ICC,KEY,NCALL)
C        THIS SUBROUTINE APPLIES A NEWTON-RAPHSON ROOT-FINDING
C        TECHNIQUE TO FIND A VALUE OF EPSILON FOR A PARTICULAR USE
C        OF THE EXPONENTIAL STRECHING TRANSFORMATION.
      DIMENSION R(40)
      FMXL=FMX
      FMINL=FMIN
      DFML=DFM
      FPCCL=FPCC
      ICCL=ICC
      GO TO (1,2),KEY
1     FNPTM2=NPT-2
      IF(NCALL.EQ.1)  EPS=(FMXL/DFML)**(1.0/FNPTM2)-1.0
      DO 3 NIT=1,ICCL
      EP1=EPS+1.0
      EP1TN=EP1**FNPTM2
      REPS=1.0/EPS
      DFMOE=DFML*REPS
      F=FMXL-FMINL-DFMOE*(EP1TN*EP1-1.0)
      IF(ABS(F).LT.FPCCL)  GO TO 4
      DFMOE2=DFMOE*REPS
      FPN=DFMOE2*(1.0+EP1TN*(EPS*FNPTM2-1.0))
      EPS=EPS+F/FPN
3     CONTINUE
      GO TO 5
2     NPTM=NPT-1
      IF(NCALL.EQ.1)
     1  EPS=((FMXL/DFML)**(1.0/(NPT-2))-1.0)*SQRT(FLOAT(NPTM))
      DO 6 L=1,NPTM
6     R(L)=1.0/SQRT(FLOAT(L))
      DO 7 NIT=1,ICCL
      SUM1=0.0
      SUM2=0.0
      DO 8 L=1,NPTM
      FLM2=L-2
      FACT1=1.0+EPS*R(L)
      FACT2=FACT1**FLM2
      SUM1=SUM1+FACT2*FACT1
8     SUM2=SUM2+(L-1)*FACT2*R(L)
      F=FMXL-FMINL-DFML*SUM1
      IF(ABS(F).LT.FPCCL)  GO TO 4
      FPN=DFML*SUM2
      EPS=EPS+F/FPN
7     CONTINUE
5     EPSIL=EPS
      WRITE(6,100)
      RETURN
4     EPSIL=EPS
      WRITE(6,101)  EPSIL,F,NIT
      RETURN
100   FORMAT(/42H EXCEEDED MAX. NO. OF ITERATIONS IN EPSIL.)
101   FORMAT(/7H EPSIL=,F12.5,5X,7H AND F=,F12.5,5X,7H AFTER ,I3,
     *  12H ITERATIONS.)
      END
*DECK FCT
      SUBROUTINE FCT(A,B,XN,IW,X,Y)
      IF(IW.EQ.1) Y=B*(1.0-(X/A)**XN)**(1.0/XN)
      IF(IW.EQ.2) X=A*(1.0-(Y/B)**XN)**(1.0/XN)
      RETURN
      END
*DECK FILLN
      SUBROUTINE FILLN(JMAX,JOMAX,KMAX,KOMAX,WB,WT,WOET,X,IW)
      DIMENSION X(40,40),A(40),B(40),C(40),D(40),F(40),WB(40),WT(40),
     1                                                          WOET(40)
      JMM=JOMAX-1
      KMM=KOMAX-1
      DO 500 K=2,KMM
      J1=2
      IF((KOMAX.GT.KMAX.AND.K.LE.KMAX+1).OR.IW.LT.0) J1=JMAX+2
      DO 100 J=J1,JMM
      WK2=.5*(WB(J)+(K-1)*(WT(J)-WB(J))/KMM)
      IF(JOMAX.GT.JMAX.AND.J.LE.JMAX+1) WK2=.5*(WOET(J)+(K-KMAX-1)*
     1                                       (WT(J)-WOET(J))/(KMM-KMAX))
      A(J)=(1.0-WK2)
      B(J)=-2.0
      C(J)=(1.0+WK2)
      F(J)=0.0
  100 CONTINUE
      F(J1)=-A(J1)*X(J1-1,K)
      IF(IABS(IW).EQ.2) F(J1)=-A(J1)*X(K,J1-1)
      F(JMM)=-C(JMM)*X(JOMAX,K)
      IF(IW.EQ.2) F(JMM)=-C(JMM)*X(K,JOMAX)
      CALL TRIB(A,B,C,D,F,J1,JMM)
      IF(IABS(IW).EQ.2) GO TO 300
      DO 200 J=J1,JMM
  200 X(J,K)=F(J)
      GO TO 500
  300 DO 400 J=J1,JMM
  400 X(K,J)=F(J)
  500 CONTINUE
      RETURN
      END
*DECK INPUT
      SUBROUTINE INPUT(X,Y,JMAX,KMAX,JOMAX,KOMAX,MAXIT,NCASE,XN,XNO,
     1IPIN,XMAX,XOMIN,XOMAX,YMAX,YOMIN,YOMAX,DX1,DXO1,DY1,DYO1,DS1,DSO1,
     2                 OMEGA,IB,IBLOUP,IWRITE,IPRINT,XVAL,ITITLE,ISTOP)
C     THE CARTESIAN COORDINATES OF THE NOZL3D CODE ARE X,Y, AND Z.
C     X IS ORIENTED IN THE GENERAL STREAMWISE DIRECTION.  IN THIS CODE-
C     RGRIDD- X AND Y ARE THE EQUIVALENT OF THE Y AND Z CARTESIAN COORDINATES
C     IN THE NOZL3D CODE.
C     JMAX,KMAX ARE # OF POINTS IN X,Y DIRECTION OF INTERIOR NOZZLE GRID
C     JOMIN,JOMAX,KOMIN,KOMAX ARE INDEX LIMITS OF OUTER GRID
C     MAXIT=MAX NUMBER OF ITERATIONS IN RELAX
C     NCASE=1(2) MEANS CALCULATE GRID FOR INSIDE(AND OUTSIDE) OF NOZZLE
C     XN(XNO) IS THE SUPER ELLIPSE NUMBER FOR THE INSIDE(OUTSIDE)
C     SURFACE OF THE NOZZLE WALL
C     IPIN=0(1) DO NOT(DO) PLOT INITIAL CONDITIONS
C     XMAX,YMAX= SEMIAXES A,B OF SUPER ELLIPSE FOR NOZZLE INTERIOR WALL
C     XOMIN,YOMIN= SEMIAXES A,B OF SUPER ELLIPSE FOR NOZZLE OUTER WALL
C     XOMAX(YOMAX)=X(Y) COORDINATE OF OUTSIDE LIMITS OF GRID ALONG X(Y) AXIS
C     XVAL= LOCATION OF THE NOZL3D X COORDINATE (FOR OUTPUT PURPOSES ONLY)
C     XVAL REPRESENTS THE NOZL3D STREAMWISE X STATION AT WHICH THE
C     CROSS SECTIONAL GRID IS GENERATED
C     DX1(DY1)=DX(DY) OF GRID NEAREST INNER NOZZLE WALL
C     DXO1(DYO1)=DX(DY) OF GRID NEAREST OUTER NOZZLE WALL
C     DS1(DSO1)=DELTA ARC LENGTH AT CORNER FOR INSIDE(OUTSIDE) SUPER ELLIPSE
C     OMEGA=OVER-RELAXATION FACTOR
C     IB(I)=0(1)[2]  DONT(DO) READ IN THE GRID POINTS ON THE ITH BOUNDRY
C                    [DISTRIBUTE PTS ON IB(3) & IB(4) TO CONFORM TO THE
C                    DISTRIBUTION ON IB(6,7) & IB(2,8)]
C     IB(I) = 4  MEANS DO ELLIPSE TEST CASE
C
C                                IB(4)
C
C  YOMAX       ---------------------------------------------
C              1                                           1
C              1                                           1
C              1                                           1
C              1                                           1
C              1                                           1
C              1                                           1
C        IB(6) 1                                           1
C              1                                           1
C              1                                           1
C              1                                           1
C              1                                           1
C              1                                           1
C              1                                           1
C              1        IB(8)                              1
C   YOMIN      1----------------------                     1 IB(3)
C    YMAX      1---------------------1                     1
C              1        IB(10)      11                     1
C              1                    11                     1
C              1                    11                     1
C              1                    11                     1
C              1                    11                     1
C              1                    11                     1
C        IB(5) 1              IB(9) 11 IB(7)               1
C              1                    11                     1
C              1                    11                     1
C              1                    11                     1
C              1                    11                     1
C              1                    11                     1
C              1                    11                     1
C              1                    11                     1
C              ---------------------------------------------
C              0                XMAX XOMIN              XOMAX
C
C                       IB(1)                IB(2)
C
C     IBLOUP=0(N)  PLOT A BLOW-UP OF N LINES AROUND THE NOZZLE CORNER
C     IWRITE=0(1) MEANS DONT(DO) WRITE OUT GRID ON UNIT 2
C     IPRINT=0(1) MEANS DONT(DO) PRINT OUT GRID
C     ITITLE=TITLE
C     ISTOP=0(1) MEANS READ(DONT READ) ANOTHER CASE
      DIMENSION ITITLE(13),X(40,40),Y(40,40),IB(10)
      READ(5,100) ITITLE
  100 FORMAT(13A6)
      WRITE(6,101) ITITLE
  101 FORMAT(1H1,13A6)
      READ(5,200) JMAX,KMAX,JOMAX,KOMAX,MAXIT,NCASE,IPIN,IBLOUP,
     1                                               IWRITE,IPRINT,ISTOP
  200 FORMAT(16I5)
      WRITE(6,201) JMAX,KMAX,JOMAX,KOMAX,MAXIT,NCASE,IPIN,IBLOUP,
     1                                               IWRITE,IPRINT,ISTOP
  201 FORMAT(// 92H     JMAX    KMAX    JOMAX  KOMAX   MAXIT    NCASE
     1  IPIN  IBLOUP  IWRITE  IPRINT  ISTOP  /16I8)
      READ(5,300) XMAX,XOMIN,XOMAX,YMAX,YOMIN,YOMAX,XVAL
  300 FORMAT(8F10.5)
      WRITE(6,301) XMAX,XOMIN,XOMAX,YMAX,YOMIN,YOMAX,XVAL
  301 FORMAT(//5X,4HXMAX,11X,5HXOMIN,10X,5HXOMAX,10X,4HYMAX
     1,11X,5HYOMIN,10X,5HYOMAX,11X,4HXVAL,/8E15.8)
      READ(5,300) DX1,DXO1,DY1,DYO1,DS1,DSO1,OMEGA
      WRITE(6,302) DX1,DXO1,DY1,DYO1,DS1,DSO1,OMEGA
  302 FORMAT(//6X,3HDX1,12X,4HDXO1,11X,3HDY1,12X,4HDY01,11X,3HDS1,12X,
     14HDSO1,12X,5HOMEGA/8E15.8)
      READ(5,300) XN,XNO
      WRITE(6,303) XN,XNO
  303 FORMAT(//6X,2HXN,13X,3HXNO/2E15.8)
      READ(5,200) IB
      WRITE(6,400) IB
  400 FORMAT(//3H IB,10I5)
      JP1=JMAX+1
      KP1=KMAX+1
      IF(IB(1).EQ.1) READ(5,300) (X(J,1),Y(J,1),J=1,JMAX)
      IF(IB(2).EQ.1) READ(5,300) (X(J,1),Y(J,1),J=JP1,JOMAX)
      IF(IB(3).EQ.1) READ(5,300) (X(JOMAX,K),Y(JOMAX,K),K=1,KOMAX)
      IF(IB(4).EQ.1) READ(5,300) (X(J,KOMAX),Y(J,KOMAX),J=1,JOMAX)
      IF(IB(5).EQ.1) READ(5,300) (X(1,K),Y(1,K),K=1,KMAX)
      IF(IB(6).EQ.1) READ(5,300) (X(1,K),Y(1,K),K=KP1,KOMAX)
      IF(IB(7).EQ.1) READ(5,300) (X(JP1,K),Y(JP1,K),K=1,KP1)
      IF(IB(8).EQ.1) READ(5,300) (X(J,KP1),Y(J,KP1),J=1,JP1)
      IF(IB(9).EQ.1) READ(5,300) (X(JMAX,K),Y(JMAX,K),K=1,KMAX)
      IF(IB(10).EQ.1)READ(5,300) (X(J,KMAX),Y(J,KMAX),J=1,JMAX)
      IF(NCASE.EQ.2) RETURN
      JOMAX=JMAX
      KOMAX=KMAX
      XOMAX=XMAX
      YOMAX=YMAX
      RETURN
      END
*DECK MAIN
      PROGRAM MAIN(INPUT,OUTPUT,TAPE5=INPUT,TAPE6=OUTPUT,TAPE2)
      DIMENSION X(40,40),Y(40,40),WL(40),WR(40),WB(40),WT(40),WOL(40),
     1       WOR(40),WOB(40),WOT(40),WOER(40),WOET(40),IB(10),ITITLE(13)
   10 CALL INPUT(X,Y,JMAX,KMAX,JOMAX,KOMAX,MAXIT,NCASE,XN,XNO,IPIN,XMAX,
     1XOMIN,XOMAX,YMAX,YOMIN,YOMAX,DX1,DXO1,DY1,DYO1,DS1,DSO1,OMEGA,IB,
     2                           IBLOUP,IWRITE,IPRINT,XVAL,ITITLE,ISTOP)
      ISW=1
      CALL BDRYS(1,JMAX,1,KMAX,0.0,XMAX,0.0,YMAX,XN,DX1,DY1,DS1,XMAX
     1                         ,YMAX,X,Y,IB,D1,D2,WL,WR,WB,WT,WOER,WOET)
      IF(NCASE.EQ.2) CALL BDRYS(JMAX+1,JOMAX,KMAX+1,KOMAX,XOMIN,XOMAX,
     1YOMIN,YOMAX,XNO,DXO1,DYO1,DSO1,XOMIN,YOMIN,X,Y,IB,D1,D2,WOL,WOR,
     2                                               WOB,WOT,WOER,WOET)
C     CALL FILL(NCASE,JMAX,JOMAX,KMAX,KOMAX,X,Y)
      CALL FILLN (JMAX,JMAX,KMAX,KMAX,WB,WT,WOET,X,1)
      CALL FILLN (KMAX,KMAX,JMAX,JMAX,WL,WR,WOER,Y,2)
      IF(NCASE.EQ.2) CALL FILLN(JMAX,JOMAX,KMAX,KOMAX,WOB,WOT,WOET,X
     1   ,1)
      IF(NCASE.EQ.2) CALL FILLN(KMAX,KOMAX,JMAX,JOMAX,WOL,WOR,
     1  WOER,Y,2)
      ITITLE(5)=6HINITIA
      ITITLE(6)=6HL COND
      ITITLE(7)=6HITIONS
      IF(IPIN.EQ.1) CALL PLAWT(JOMAX,KOMAX,X,Y,XOMAX,0.0,YOMAX,0.0
     1      ,XNO,DY1,DX1,DYO1,DXO1,DS1,ITITLE,JMAX,KMAX,XN,JOMAX,KOMAX)
      CALL RELAX(JMAX,KMAX,JMAX,KMAX,X,Y,OMEGA,MAXIT,WB,WT,WL,WR,WOER,
     1                                                            WOET)
C     CALL RELAXS(JMAX,KMAX,X,Y,OMEGA,MAXIT,WL,WR,WB,WT)
      IT4=4*MAXIT
      IF(ISW.EQ.-1) IT4=MAXIT
      IF(NCASE.EQ.2) CALL RELAX(JMAX,KMAX,JOMAX,KOMAX,X,Y,OMEGA,IT4,
     1                                        WOB,WOT,WOL,WOR,WOER,WOET)
      ITITLE(5)=6H FINAL
      ITITLE(6)=6H GRID
      ITITLE(7)=6H
      CALL PLAWT(JOMAX,KOMAX,X,Y,XOMAX,0.,YOMAX,0.
     1    ,XNO,DY1,DX1,DYO1,DXO1,DS1,ITITLE,JMAX,KMAX,XN,JOMAX,KOMAX)
      IF(IWRITE.EQ.1) WRITE(2)XVAL,((X(J,K),Y(J,K),J=1,JOMAX),K=1,KOMAX)
      IF(IPRINT.EQ.1) CALL OUTPUT(X,Y,JOMAX,KOMAX,XVAL)
C     IF(NCASE.EQ.2) GO TO 100
C     CALL CLUSTR(JMAX,X,Y,D1,D2,KMAX)
C     ITITLE(5)=6HEXPONE
C     ITITLE(6)=6HNTIALL
C     ITITLE(7)=6HY CLUS
C     ITITLE(8)=6HTERRED
C     ITITLE(9)=6H GRID
C     CALL PLAWT(JOMAX,KOMAX,X,Y,XOMAX,0.0,YOMAX,0.0
C    1       ,XNO,DY1,DX1,DYO1,DXO1,DS1,ITITLE,JMAX,KMAX,XN,JOMAX,KOMAX)
C 100 CONTINUE
      IF(IBLOUP.GT.0)CALL PLOTBU(JOMAX,KOMAX,X,Y,XNO,XN,DY1,DX1,DYO1,DXO1
     1                                     ,DS1,ITITLE,JMAX,KMAX,IBLOUP)
      IF(ISTOP.EQ.0) GO TO 10
      CALL EOFTV
      STOP
      END
*DECK OUTER
      SUBROUTINE OUTER(XMAX,XMIN,YMAX,YMIN,XORG,ETAC,BETA,JMAX,KMAX,X,Y)
C
C        THIS SUBROUTINE PLACES POINTS ON BOTTOM-FRONT-TOP BOUNDARY IN
C        ANGULAR FASHION.
C
      DIMENSION X(40,40),Y(40,40)
C
      LOGICAL CLUSTR
C
      DATA PI/3.141592654/
C
      SINH(X)=0.5*(EXP(X)-EXP(-X))
C
      ETARU=ATAN2(YMAX,XMAX-XORG)
      ETARL=ATAN2(-YMIN,XMAX-XORG)
      DETA=(2.0*PI-(ETARU+ETARL))/(JMAX-1)
      CLUSTR = .FALSE.
      IF( BETA.GT.0.0) CLUSTR = .TRUE.
      IF(.NOT. CLUSTR) GO TO 14
      FACT=PI/(2.0*PI-(ETARU+ETARL))
      FACTR=1.0/FACT
      ETACT=(ETAC-ETARU)*FACT
      B=0.5*ALOG((1.+(EXP(BETA)-1.)*ETACT/PI)/(1.+(EXP(-BETA)-1.)*
     1  ETACT/PI))
      RSB = 1./SINH(B)
14    ETA=ETARU
      ANG1=ATAN2(YMAX,XMIN-XORG)
      ANG2=ATAN2(YMIN,XMIN-XORG)+2.0*PI
      NSIDE = 1
C
      ETARUD=ETARU*180./PI
      ETARLD=ETARL*180./PI
      ANG1D=ANG1*180./PI
      ANG2D=ANG2*180./PI
      WRITE(6,109)  CLUSTR,ETARUD,ETARLD,ANG1D,ANG2D
      IF(CLUSTR.AND.ETAC.LT.ETARU.OR.CLUSTR.AND.ETAC.GT.(2.0*PI-ETARL))
     1  GO TO 22
C
      DO 9 JJ=2,JMAX
      J=JMAX+1-JJ
      ETA = ETA + DETA
      IF(.NOT.CLUSTR)  GO TO 26
      ETAT=(ETA-ETARU)*FACT
      PHIT=ETACT*(SINH(BETA*ETAT/PI-B)*RSB+1.)
      PHI=ETARU+PHIT*FACTR
      GO TO 27
26    PHI=ETA
C
27    GO TO (1,2,3),NSIDE
1     IF(PHI.GT.ANG1)  NSIDE=2
      GO TO 3
2     IF(PHI.GT.ANG2)  NSIDE=3
C
3     GO TO (10,11,12),NSIDE
   10 Y(J,KMAX) = YMAX
      X(J,KMAX) = YMAX/TAN(PHI) + XORG
      GO TO 21
   11 X(J,KMAX) = XMIN
      Y(J,KMAX) = ( XMIN - XORG)*TAN(PHI)
      GO TO 21
12    X(J,KMAX)=XORG+YMIN/TAN(PHI)
      Y(J,KMAX)=YMIN
21    CONTINUE
C
      ETAD=ETA*180./PI
      PHID=PHI*180./PI
      WRITE(6,113)  J,KMAX,X(J,KMAX),J,KMAX,Y(J,KMAX),ETAD,PHID,NSIDE
    9 CONTINUE
C
      RETURN
C
22    WRITE(6,114)
      STOP
C
109   FORMAT(/8H CLUSTR=,L1,5X,6HETARU=,F8.2,5X,6HETARL=,F8.2,5X,
     1  5HANG1=,F8.2,5X,5HANG2=,F8.2//24H OUTER BOUNDARY ON TOP, ,
     2 18HFRONT, AND BOTTOM&)
113   FORMAT(3H X(,I3,1H,,I3,2H)=,F12.5,5X,2HY(,I3,1H,,I3,2H)=,F12.5,
     1  5X,4HETA=,F8.2,5X,4HPHI=,F8.2,5X,6HNSIDE=,I1)
114   FORMAT(/23H ERROR EXIT.  BAD ETAC.)
C
      END
*DECK OUTPUT
      SUBROUTINE OUTPUT(X,Y,JMAX,KMAX,XVAL)
      DIMENSION X(40,40),Y(40,40)
      WRITE(6,100) XVAL
  100 FORMAT(1H1,5X,3HX =,F10.5///5X,7HY ARRAY/)
      IF(JMAX.LE.20) GO TO 400
      WRITE(6,200) ((X(J,K),J=1,JMAX),K=1,KMAX)
  200 FORMAT(20F6.3/3X,20F6.3)
      WRITE(6,300)
  300 FORMAT(///5X,7HZ ARRAY)
      WRITE(6,200) ((Y(J,K),J=1,JMAX),K=1,KMAX)
      RETURN
  400 WRITE(6,500) ((X(J,K),J=1,JMAX),K=1,KMAX)
  500 FORMAT(20F6.3)
      WRITE(6,300)
      WRITE(6,500) ((Y(J,K),J=1,JMAX),K=1,KMAX)
      RETURN
      END
*DECK PLAWT
      SUBROUTINE PLAWT(N,M,X,Y,XMAX,XMIN,YMAX,YMIN,XNN,DY1,DX1,
     1            DYO1,DXO1,DS1,ITITLE,JMAX,KMAX,XNI,JOMAX,KOMAX)
C
C        THIS SUBROUTINE IS ONE OF THE SUBROUTINES THAT PLOTS GRIDS
C        ON THE SC4020 PLOTTER.
C
      COMMON/PLOTC/SING
      DIMENSION X(40,40),Y(40,40),XX(40),YY(40),TIT(15),ITITLE(13)
C
C        READJUST PLOT LIMITS SO AS TO AVIOD A STRECHED PLOT.
      IW=1
      IF(JOMAX.GT.JMAX.AND.KOMAX.EQ.KMAX) IW=0
      XDIF=XMAX-XMIN
      YDIF=YMAX-YMIN
      IF(XDIF.LT.YDIF)  GO TO 4
      XDIFH=XDIF*0.5
      YMID=(YMAX+YMIN)*0.5
      YMX=YMID+XDIFH
      YMN=YMID-XDIFH
      XMX=XMAX
      XMN=XMIN
      GO TO 5
4     YDIFH=YDIF*0.5
      XMID=(XMAX+XMIN)*0.5
      XMX=XMID+YDIFH
      XMN=XMID-YDIFH
      YMX=YMAX
      YMN=YMIN
5     CONTINUE
C
C        PLOT THE LINES.
      CALL AXIS(XMN,XMX,YMN,YMX,0,0,XMIN,XMAX,YMIN,YMAX)
      CALL TITLE(ITITLE,78)
      ENCODE(63,10,TIT) M,N,XNN,DY1,DX1,DS1
   10 FORMAT(5HKMAX=,I2,6H JMAX=,I2,4H XN=,F5.1,5H DY1=,F8.7,5H DX1=
     1                                               ,F8.7,5H DS1=,F8.7)
      CALL PRINTV(63,TIT,80,950)
      IF(M.EQ.KMAX) GO TO 30
      ENCODE(27,20,TIT) KMAX,JMAX,XNI
   20 FORMAT(6HKIMAX=,I2,7H JIMAX=,I2,5H XNI=,F5.1)
      CALL PRINTV(27,TIT,150,925)
   30 DO1J=1,M
      MAX=JMAX
      IF(J.GT.KMAX.OR.IW.EQ.0) MAX=N
1     CALL PLOT(X(1,J),Y(1,J),N,0,MAX)
      DO2I=1,N
      DO3J=1,M
      XX(J)=X(I,J)
    3 YY(J)=Y(I,J)
      MAX=KMAX
      IF(I.GT.JMAX) MAX=M
2     CALL PLOT(XX,YY,M,0,MAX)
      IF(SING.LE.0.0) RETURN
      IX1=IXV(X(JMAX,1))
      IY1=IYV(Y(JMAX,1))
      CALL LINEV(IX1,IY1-15,IX1,IY1)
      CALL LINEV(IX1,IY1-15,IX1,IY1)
      RETURN
      END
*DECK PLOT
      SUBROUTINE PLOT(X,Y,NBR,NSYM,MAX)
C
C        THIS SUBROUTINE IS ONE OF THOSE THAT PLOTS GRIDS
C        ON THE SC4020 PLOTTER.
C
      DIMENSION X(1),Y(1),MARKPT(5)
      DATA MARKPT /42,16,55,38,44/
C     SYMBOLS ARE, IN ORDER&  . + X O *
      IF (NSYM .GT. 0 .AND. NSYM .LE. 5) GO TO 100
      J=IABS(NBR)-1
      DO 110 I=1,J
      IF(I.EQ.MAX) GO TO 110
           CALL LINEV( IXV(X(I)), IYV(Y(I)), IXV(X(I+1)),IYV(Y(I+1)) )
 110  CONTINUE
      RETURN
 100  CALL APLOTV(NBR,X,Y,1,1,1,MARKPT(NSYM))
      RETURN
      END
*DECK PLOTBU
      SUBROUTINE PLOTBU(JOMAX,KOMAX,X,Y,XNO,XN,DY1,DX1,DYO1,DXO1,DS1,
     1                                          ITITLE,JMAX,KMAX,IBLOUP)
      DIMENSION X(40,40),Y(40,40),ITITLE(1)
      IF(JMAX.EQ.JOMAX) GO TO 100
      IBH=IBLOUP/2
      J1=JMAX-IBH+1
      J2=JMAX+IBH
      K1=KMAX-IBH+1
      K2=KMAX+IBH
      MAXJ=IBH
      MAXK=IBH
      GO TO 200
  100 J1=MAX0(JMAX-2*IBH+1,1)
      J2=JMAX
      K1=MAX0(KMAX-2*IBH+1,1)
      K2=KMAX
      MAXJ=JMAX
      MAXK=KMAX
  200 JJ=0
      XMIN=1.0E36
      YMIN=1.0E36
      XMAX=-1.0E36
      YMAX=-1.0E36
      DO 300 J=J1,J2
      JJ=JJ+1
      KK=0
      DO 300 K=K1,K2
      KK=KK+1
      X(JJ,KK)=X(J,K)
      Y(JJ,KK)=Y(J,K)
      XMIN=AMIN1(XMIN,X(JJ,KK))
      XMAX=AMAX1(XMAX,X(JJ,KK))
      YMIN=AMIN1(YMIN,Y(JJ,KK))
      YMAX=AMAX1(YMAX,Y(JJ,KK))
  300 CONTINUE
      ITITLE(5)=6HBLOWUP
      ITITLE(6)=6H OF MO
      ITITLE(7)=6HST CON
      ITITLE(8)=6HGESTED
      ITITLE(9)=6H AREA
      CALL PLAWT(JJ,KK,X,Y,XMAX,XMIN,YMAX,YMIN,XNO,DY1,DX1,
     1                 DY01,DX01,DS1,ITITLE,MAXJ,MAXK,XN,JOMAX,KOMAX)
      RETURN
      END
*DECK RELAX
      SUBROUTINE RELAX(JMAX,KMAX,JOMAX,KOMAX,X,Y,OMEGA,MAXIT,WB,WT,WL,WR
     1                                                       ,WOER,WOET)
C        THIS SUBROUTINE SOLVES BY SLOR THE DIFFERENTIAL EQUATIONS BASED
C        ON A MODIFIED THOMPSON-THAMES-MASTIN S METHOD OF GENERATING GRIDS.
C     THE MODIFICATION EMPLOYS SPECIALIZED SOURCE TERMS INVOLVING FREE
C     PARAMETERS WP AND WQ THAT ARE COMPUTED FROM THE BOUNDARY VALUES.
      DIMENSION X(40,1),WB(1),WT(1),WL(1),WR(1),WOER(1),WOET(1),Y(40,1)
      DIMENSION A(40),B(40),C(40),D(40),F(40),G(40)
      IW=1
      IF(JOMAX.GT.JMAX.AND.KOMAX.EQ.KMAX) IW=0
      KMM=KOMAX-1
      JMM=JOMAX-1
      ICOUNT=0
2     ICOUNT=ICOUNT+1
      RSUM=0.
      RXSUM=0.0
      RYSUM=0.0
      DO 1 K=2,KMM
      J1=2
      IF((KOMAX.GT.KMAX.AND.K.LE.KMAX+1).OR.IW.EQ.0) J1=JMAX+2
      DO 3 J=J1,JMM
      XXD=(X(J+1,K)-X(J-1,K))*0.5
      XED=(X(J,K+1)-X(J,K-1))*0.5
      YXD=(Y(J+1,K)-Y(J-1,K))*0.5
      YED=(Y(J,K+1)-Y(J,K-1))*0.5
      AD=XED**2+YED**2
      BD=XXD*XED+YXD*YED
      GD=XXD**2+YXD**2
      XXED=(X(J+1,K+1)-X(J+1,K-1)-X(J-1,K+1)+X(J-1,K-1))*0.25
      YXED=(Y(J+1,K+1)-Y(J+1,K-1)-Y(J-1,K+1)+Y(J-1,K-1))*0.25
      BD=-2.0*BD
      WQ=.5*(WL(K)+(J-1.0)*(WR(K)-WL(K))/JMM)
      WP=.5*(WB(J)+(K-1.0)*(WT(J)-WB(J))/KMM)
      IF(JOMAX.EQ.JMAX.OR.IW.EQ.0) GO TO 20
      IF(J.LE.JMAX+1)WP=.5*(WOET(J)+(K-KMAX-1)*(WT(J)-WOET(J))/(KMM-KMAX
     1                                                                ))
      IF(K.LE.KMAX+1)WQ=.5*(WOER(K)+(J-JMAX-1)*(WR(K)-WOER(K))/(JMM-JMAX
     1                                                                ))
   20 A(J)=AD*(1.0-WP)
      B(J)=-AD-AD-GD-GD
      C(J)=AD*(1.0+WP)
      F(J)=-BD*XXED-GD*((1.0+WQ)*X(J,K+1)+(1.0-WQ)*X(J,K-1))
3     G(J)=-BD*YXED-GD*((1.0+WQ)*Y(J,K+1)+(1.0-WQ)*Y(J,K-1))
      F(J1)=F(J1)-A(J1)*X(J1-1,K)
      G(J1)=G(J1)-A(J1)*Y(J1-1,K)
      F(JMM)=F(JMM)-C(JMM)*X(JOMAX,K)
      G(JMM)=G(JMM)-C(JMM)*Y(JOMAX,K)
      CALL TRIB(A,B,C,D,F,J1,JMM)
      CALL TRIB(A,B,C,D,G,J1,JMM)
      DO 4 J=J1,JMM
      XC=OMEGA*(F(J)-X(J,K))
      YC=OMEGA*(G(J)-Y(J,K))
      RSUM=RSUM+ABS(XC)+ABS(YC)
      RXSUM=RXSUM+ABS(XC)
      RYSUM=RYSUM+ABS(YC)
      X(J,K)=X(J,K)+XC
4     Y(J,K)=Y(J,K)+YC
    1 CONTINUE
      WRITE(6,100)RSUM,RXSUM,RYSUM,ICOUNT
      IF(ICOUNT.LT.MAXIT)  GO TO  2
      RETURN
100   FORMAT(29H SUM OF RESIDUALS (X+Y),X,Y = ,3F20.10,
     1 7H AFTER ,I5,12H ITERATIONS.)
      END
*DECK ROB
      SUBROUTINE ROB(X,Y,WR,JMIN,JMAX,KMIN,KMAX,YMAX)
      DIMENSION X(40,40),Y(40,40)
      YM=.5*(Y(1,KMIN+1)+Y(JMIN,KMIN))
      DYMB=Y(JMIN,KMIN)-Y(JMIN,KMIN-1)
      DYMT=Y(1,KMIN+1)-Y(1,KMIN)
      DYM=.5*(DYMB+DYMT)
      Y(JMAX,KMIN)=YM
      X(JMAX,KMIN)=X(JMAX,1)
      Y(JMAX,KMIN-1)=YM-DYM
      X(JMAX,KMIN-1)=X(JMAX,1)
      Y(JMAX,KMIN+1)=YM+DYM
      X(JMAX,KMIN+1)=X(JMAX,1)
      Y(JMAX,KMAX)=YMAX
      X(JMAX,KMAX)=X(JMAX,1)
      DYM=1.1*DYM
      E=EPSIL(Y(JMAX,KMAX)-Y(JMAX,KMIN+1),0.,DYM,KMAX-KMIN,.001,100,1,1)
      K1=KMIN+2
      KMAXM1=KMAX-1
      DO 100 K=K1,KMAXM1
      X(JMAX,K)=X(JMAX,1)
  100 Y(JMAX,K)=Y(JMAX,K-1)+DYM*(1.0+E)**(K-K1)
      E=EPSIL(Y(JMAX,KMIN-1)-Y(JMAX,1),0.0,DYM,KMIN-1,.001,100,1,1)
      K2=KMIN-2
      DO 200 KK=2,K2
      K=K2-KK+2
      X(JMAX,K)=X(JMAX,1)
  200 Y(JMAX,K)=Y(JMAX,K+1)-DYM*(1.0+E)**(K2-K)
      CALL WCALC(X,Y,WR,JMAX,JMAX,2,KMAXM1)
      RETURN
      END
*DECK SLOPE
      SUBROUTINE SLOPE(A,B,XN,X2,Y2)
      JS=1000
      DX=A/(JS-1.0)
      X1=0.0
      Y1=B
      DO 100 J=2,JS
      X2=(J-1)*DX
      CALL FCT(A,B,XN,1,X2,Y2)
      S=(Y2-Y1)/(X2-X1)
      IF(S.LT.-1.0) RETURN
      X1=X2
      Y1=Y2
  100 CONTINUE
      WRITE(6,200)
  200 FORMAT(5X,27HSTOPPED IN SUBROUTINE SLOPE )
      STOP
      END
*DECK TITLE
      SUBROUTINE TITLE(ITITLE,NCHARS)
C
C        THIS SUBROUTINE IS ONE OF THOSE THAT PLOTS GRIDS
C        ON THE SC4020 PLOTTER.
C
C     THE MAXIMUM NUMBER OF CHARACTERS ALLOWED IN THE TITLE IS 108.
C
      ICHARS=IABS(NCHARS)
      IF (ICHARS .GT. 108) ICHARS=108
      IF (ICHARS .GT. 54) IX=14
      IF (ICHARS .LE. 54) IX=510-(ICHARS/2)*18
      IY=990
      CALL RITE2V(IX,IY,1010,90,1,ICHARS,1,ITITLE,NLAST)
      RETURN
      END
*DECK TOB
      SUBROUTINE TOB(X,Y,WT,JMIN,JMAX,KMIN,KMAX)
      DIMENSION X(40,40),Y(40,40)
      XM=.5*(X(JMIN+1,1)+X(JMIN,KMIN))
      DXMB=X(JMIN,KMIN)-X(JMIN-1,KMIN)
      DXMT=X(JMIN+1,1)-X(JMIN,1)
      DXM=.5*(DXMB+DXMT)
      X(JMIN,KMAX)=XM
      Y(JMIN,KMAX)=Y(1,KMAX)
      X(JMIN-1,KMAX)=XM-DXM
      Y(JMIN-1,KMAX)=Y(1,KMAX)
      X(JMIN+1,KMAX)=XM+DXM
      Y(JMIN+1,KMAX)=Y(1,KMAX)
      DXM=1.1*DXM
      E=EPSIL(X(JMAX,KMAX)-X(JMIN+1,KMAX),0.,DXM,JMAX-JMIN,.001,100,1,1)
      J1=JMIN+2
      JMAXM1=JMAX-1
      DO 100 J=J1,JMAXM1
      X(J,KMAX)=X(J-1,KMAX)+DXM*(1.0+E)**(J-J1)
  100 Y(J,KMAX)=Y(1,KMAX)
      E=EPSIL(X(JMIN-1,KMAX)-X(1,KMAX),0.0,DXM,JMIN-1,.001,100,1,1)
      J2=JMIN-2
      DO 200 JJ=2,J2
      J=J2-JJ+2
      X(J,KMAX)=X(J+1,KMAX)-DXM*(1.0+E)**(J2-J)
  200 Y(J,KMAX)=Y(1,KMAX)
      CALL WCALC(X,Y,WT,2,JMAX-1,KMAX,KMAX)
      RETURN
      END
*DECK TRIB
      SUBROUTINE TRIB(A,B,C,X,F,NL,NU)
      DIMENSION A(2),B(2),C(2),X(2),F(2)
C
C        THIS SUBROUTINE SOLVES A TRI-DIAGONAL SYSTEM OF LINEAR
C        EQUATIONS.
C
      X(NL)=C(NL)/B(NL)
      F(NL)=F(NL)/B(NL)
      NLP1 = NL +1
      DO 1 J=NLP1,NU
      Z=1./(B(J)-A(J)*X(J-1))
      X(J)=C(J)*Z
 1    F(J)=(F(J)-A(J)*F(J-1))*Z
      NUPNL=NU+NL
      DO 2 J1=NLP1,NU
      J=NUPNL-J1
 2    F(J)=F(J)-X(J)*F(J+1)
      RETURN
      END
*DECK WCALC
      SUBROUTINE WCALC(X,Y,W,J1,J2,K1,K2)
      DIMENSION X(40,40),Y(40,40),W(1)
      IF(J1.EQ.J2) GO TO 200
      DO 100 J=J1,J2
      T=X(J+1,K1)-X(J-1,K1)
      T2=Y(J+1,K1)-Y(J-1,K1)
      W(J)=-2.0*(T*(X(J+1,K1)-2.0*X(J,K1)+X(J-1,K1))+
     1          T2*(Y(J+1,K1)-2.0*Y(J,K1)+Y(J-1,K1)))/(T**2+T2**2)
  100 CONTINUE
      RETURN
  200 DO 300 K=K1,K2
      T=Y(J1,K+1)-Y(J1,K-1)
      T2=X(J1,K+1)-X(J1,K-1)
      W(K)=-2.0*(T*(Y(J1,K+1)-2.0*Y(J1,K)+Y(J1,K-1))+
     1          T2*(X(J1,K+1)-2.0*X(J1,K)+X(J1,K-1)))/(T**2+T2**2)
  300 CONTINUE
      RETURN
      END
*DECK XOFS
      SUBROUTINE XOFS(DS,X1,DX,A,B,XN,EPSLON)
      IF(X1+DX.GE.A) DX=.999999*(A-X1)
      H=-.1*DX
      CALL ARCLEN(X1,X1+DX,A,B,XN,16,DSDX)
      DXP=0.0
      DXM=A
      IF(DS-DSDX.GT.0.0) DXP=DX
      IF(DS-DSDX.LT.0.0) DXM=DX
      DO 500 I=1,20
      IF(X1+DX+H.GE.A) H=-.9*(A-X1-DX)
      IF(DX+H.LT.0.0) H=-DX/2.0
      CALL ARCLEN(X1,X1+DX+H,A,B,XN,16,DSDXPH)
      D=(DSDXPH-DSDX)/H
      DX=DX-(DSDX-DS)/D
      IF(X1+DX.GE.A) DX=.999999*(A-X1)
      IF(DXP.EQ.0.0.OR.DXM.EQ.A) GO TO 450
      IF(DX.GE.DXP.OR.DX.LE.DXM) DX=(DXP+DXM)/2.0
  450 CALL ARCLEN(X1,X1+DX,A,B,XN,16,DSDX)
      IF(ABS(DS-DSDX)/DS.LT.EPSLON) GO TO 600
      IF(DS-DSDX.GT.0.0) DXP=AMAX1(DXP,DX)
      IF(DS-DSDX.LT.0.0) DXM=AMIN1(DXM,DX)
      H=(DSDX-DS)/D
  500 CONTINUE
      WRITE(6,550) X1,DX,H,D,DSDX,DSDXPH,DS,DXP,DXM,A
  550 FORMAT(5X,26HSTOPPED IN SUBROUTINE XOFS /10E13.7)
      STOP
  600 X1=X1+DX
      RETURN
      END
